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Abstract 



In this article, we analyse three related preconditioned steepest descent algo- 
rithms, which are partially popular in Hartree-Fock and Kohn-Sham theory as well 
as invariant subspace computations, from the viewpoint of minimization of the cor- 
responding functionals, constrained by orthogonality conditions. We exploit the ge- 
ometry of the of the admissible manifold, i.e. the invariancc with respect to unitary 
^ transformations, to reformulate the problem on the Grassmann manifold as the ad- 

' missible set. We then prove asymptotical linear convergence of the algorithms under 

^ the condition that the Hessian of the corresponding Lagrangian is elliptic on the 

tangent space of the Grassmann manifold at the minimizer. 

as 

in 
o 

g 1 Introduction 

> 
>< 



On the length-scale of atomistic or molecular systems, physics is governed by the laws of 
quantum mechanics. A reliable computation required in various fields in modern sciences 
and technology should therefore be based on the first principles of quantum mechanics, so 
that ab initio computation of the electronic wave function from the stationary electronic 
Schrodinger equation is a major working horse for many applications in this area. To 
reduce computational demands, the high dimensional problem of computing the wave 
function for electrons is often, for example in Hartree-Fock and Kohn-Sham theory, 
replaced by a nonlinear system of equations for a set $ = {(pi, . . . ^ip^) of single particle 
wave functions 9?i(x) eV = H^{E?). This ansatz corresponds to the following abstract 
formulation for the minimization of a suitable energy functional J: 



*This work was supported by the DFG SPP 1445: "Modern and universal first-principles methods for 
many-electron systems in chemistry and physics" and the EU NEST project BigDFT. 
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Problem 1: Minimize 



J-.V^'^R, J{<i>) = J{^i,...,^n) ^ min, (1.1) 
where JT" is a sufficiently often different iable functional which is 

(i) invariant with respect to unitary transformations, i.e. 

N 

= j($u) = (1.2) 

for any orthogonal matrix U G M"^", and 

(ii) subordinated to the orthogonality constraints 

{(pi,(Pj) := / (pi{x)ipj{x)dx = 6ij. (1.3) 



R3 

In the present article, we shall be concerned with minimization techniques for J' along 



the admissible manifold characterized by (1.3). The first step towards this will be to set 



up the theoretical framework of the Grassmann manifold to be introduced in section |2j 
reflecting the constraints (i) and (ii) imposed on the functional J and the minimizer $, 
respectively. In applications in electronic structure theory, formulation of the ffist order 



optimality (necessary) condition for the problem (1.1) results in a nonlinear eigenvalue 
problem of the kind: 

A^iPi = Xi(pi, Ai < A2 < . . . < Aat (1.4) 

for N eigenvalues Aj and the corresponding solution functions assembled in $. In these 
equations, the operator A^, is a symmetric bounded linear mapping : V = H^iM.^) 
V = H~^{M.^) depending on $, so that we are in fact faced with a nonlinear eigenvalue 
problem. is called the Fock operator in Hartree-Fock theory, and Kohn-Sham Hamil- 
tonian in density functional theory (DFT) respectively. We will illustrate the relation 



between (1.4) and the minimization task above in further detail in section 3 In this work 



our emphasis will rather be on the algorithmic approximation of the minimizer of , i.e. 



an invariant subspace span[$] := span{(y9i, . . . , (fw}, of (1.4), in the corresponding energy 
space than on computation of the eigenvalues Ai, . . . , Aat. 

One possible procedure for computing the minimum of is the so-called direct mini- 
mization, utilized e.g. in DFT calculation, which performs a steepest descent algorithm 
by updating the gradient of JT, i.e. the Kohn-Sham Hamiltonian or Fock operator, in each 
iteration step. Direct minimization, as proposed in p], is prominent in DFT calculations 
if good preconditioners are available and the systems under consideration are large, e.g. 
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for the computation of electronic structure in bulk crystals using plane waves, finite dif- 
ferences [7] and the recent wavelet code developed in the BigDFT project (see |15]). In 
contrast to the direct minimization procedure is the self consistent field iteration (SCF), 
which keeps the Fock operator fixed until convergence of the corresponding eigenfunctions 
and updates the Fock operator thereafter, see section [3j 

In the rest of this article, we will pursue different variants of projected gradient algorithms 
to be compiled in section |4| In addition, we will (for the case where the gradient i7'($) 
can be written as an operator A<j, applied to $, as it is the case in electronic structure 
calculation) investigate an algorithm based on j4j following a preconditioned steepest 
descent along geodesies on the manifold, so that no re-projections onto the admissible 
manifold are required. It turns out that all these algorithms to be proposed perform in 
a similar way. For matters of rigorous mathematical analysis, let us note at this point 
that the mathematical theory about Hartree-Fock is still too incomplete to prove the 
assumptions required in the present paper; even less is known for Kohn-Sham equations, 
due to the fact that there are so many different models used in practice. If the assumptions 
are not met for a particular problem, it is not clear whether it is a deficiency of the problem 



or a real pathological situation. Along with (1.1), we will therefore consider the following 



simplified prototype problem for a fixed operator A: 
Simplified problem 2: Minimize 

N 

J7a(v5i, • • • ,V57v) := ^(V5i,^V5i) — ' min, V^j) = (1.5) 



Analogous treatment with Lagrange techniques shows that this special case of problem 1 
is the problem of computing the first eigenfunctions, resp. the lowest A^ eigenvalues of 
A (see Lemma |3]). While this is an interesting problem by itself, e.g. if A is an eigenvalue of 
multiplicity A^, it is also of interest as a sort of prototype: Properties that can be proven 
for this problem may hold in the more general case for Hartree-Fock or Kohn-Sham. In 
particular, we will show that for A symmetric and bounded from below, the Hessian of 
the Lagrangian, taken at the solution \E', is elliptic on a specific tangent manifold at \E', 
an essential ingredient to prove linear convergence of all of the proposed algorithms in 



section 5 The same convergence results will be shown to hold for (1.1) if we impose this 



ellipticity condition on the Lagrangian of of the nonlinear problem. Note that the prob- 



lem type (1.5) also arises in many other circumstances, which we will not consider here in 
detail. Let us just note that the algorithms presented in section |4] also provide reasonable 
routines for the inner cycles of the SCF procedure. 

In the context of eigenvalue computations, variants of our basic algorithm 1, applied to 
problem 2, have been considered by several authors (see e.g. [SII2H1E1]) reporting excellent 
performance, in particular if subspace acceleration techniques are applied and the precon- 
ditioner is chosen appropriately; in [iQl [TT], an adaptive variant was recently proposed 
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and analysed for the simpler case = 1. In contrast to all these papers, we will view the 
algorithms as steepest descent algorithms for optimization of J under the orthogonality 
constraints given above, as such a systematic treatment does not only simplify the proofs 
but also provides the insight necessary to understand the direct minimization techniques 
for the more complicated nonlinear problems of the kind (1.1) in DFT and HF. 
Our analysis will cover closed (usually finite dimensional) subspaces of V/j C as well 
as the energy space V itself, so that finite dimensional approximations by Ritz-Galerkin 
methods and also finite difference approximations are included in our analysis. In partic- 
ular, our results are also valid if Gaussian type basis functions are used. The convergence 
rates will be independent of the discretization parameters like mesh size. However, the 
choice of an appropriate preconditioning mapping to be used in our algorithms is crucial. 
Fortunately, such preconditioners can often easily be constructed, e.g. by the use of multi- 
grid methods for finite elements, finite differences or wavelets, polynomials [23 El C] • Our 
analysis will show that for the gradient algorithms under consideration, it suffices to use 
a fixed preconditioner respectively relaxation parameter. In particular, no expensive line 
search is required. 

All results proven will be local in nature meaning that the initial guess is supposed to 
be already sufficiently close to the exact one. At the present stage, we will for the sake 
of simplicity consider only real valued solutions for the minimization problem. Neverthe- 
less, complex valued functions can be treated by minor modifications. Note that since 
the present approach is completely based on a variational framework, i.e. considering a 
constrained optimization problem, it does not include unsymmetric eigenvalue problems 
or the computation of other eigenvalues than the lowest ones. 



2 Optimization on Grassmann manifolds 

The invariance of the functional J with respect to uniform transformations among the 
eigenfunctions shows a certain redundance inherent in the formulation of the minimization 



task (1.1). Therefore, it will be more advantageous to factor out the unitary invariance of 
the functional J7, resulting in the usage of the Stiefel and Grassmann manifolds, originally 
defined in finite dimensional Euclidean Hilbert spaces in [3], see also for an extensive 
exposition. In this section, we will generalize this concept for the present infinite dimen- 
sional space equipped with the L2 inner product. In the next section, we will then 
apply this framework to the minimization problems for the HF and KS functionals. First 
of all, we shall briefly introduce the spaces under consideration and some notations. 
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2.1 Basic notations 



Letting H = L2 := L2(M^) or a closed subspace of L2, we will work with a Gelfand triple 

V C H C V with the usual L2 inner product (.,..), as dual pairing on V x V, where either 

V := = H^{M.^) or an appropriate subspace corresponding to a Galerkin discretization. 
Because the ground state is determined by a set $ of one-particle functions ipi G V, 
we will formulate the optimization problem on an admissible subset of V^^. To this end, 
we extend inner products and operators from V to V"^ by the following 

Definitions 1. For ^ = (V'l, . . . , G 1/^, $ = {^u • • • , ^tv) e (V^)' = {VY , and the 

L2 inner product {., ..) given on H = L2, we denote 

{^^^) :=((^.,^,))f,.,GM^x^, 
and introduce the dual pairing 

N 

(($,vi/)):=tr($^vI/) = ^(^^,^^) 

1=1 

on {VY X V^. 

Because there holds = V ^ M^, we can canonically expand any operator R : V ^ V 
to an operator 

n := I : = V ^R^ V'^ , $ 7^$ = {Ripi, . . . , Rcpn). (2.1) 

Throughout this paper, for an operator V V denoted by a capital letter as A, B,D, . . ., 
the same calligraphic letter A,B,V, . . ., will denote this expansion to . 



Further, we will make use of the following operations: 

Definitions 2. For ^ e and M e M^^^, we define the set $M = (/ O M)$ 
hy ($M)j := Xlili "^«,jV^«? c/. also the notation in (1.2), and for (p E V and v = 
{vi,...,viy) G the element (j) ® v E hy (vicf), . . . ,V]\f(t>). Finally, we denote by 
0{N) the orthogonal group ofR^^^. 



2.2 Geometry of Stiefel and Grassmann manifolds 

Let us now introduce the admissible manifold and prove some of its basic properties. Note 
in this context that well established results of |1] for the case in the finite dimensional 
Euclidean spaces cannot be applied to our setting without further difficulties, because the 
norm induced by the L2 inner product is weaker than the present y-norm. 
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Our aim is to minimize the functionals where J is either Juf^ Jks or Ja-, under 

the orthogonahty constraint {ipi,ipj) = 6ij, i.e. 

= I e M^^^. (2.2) 



The subset of V satisfying the property (2.2) is called the Stiefel manifold (cf. [4 



Vv,N := {$ = e V, ($^<l>> - I = G M^^^} 

i.e. the set of all orthonormal bases of A^- dimensional subspaces of V. 



All functionals under consideration are unitarily invariant, i.e. there holds (1.2). To get 
rid of this nonuniqueness, we will identify all orthonormal bases $ G Vv,n spanning the 
same subspace V$ := span {yj, : z = 1, . . . , A^}. To this end we consider the Grassmann 
manifold, defined as the quotient 



'V,N 



of the Stiefel manifold with respect to the equivalence relation $~$ if $ = $U for any 
U G 0{N). We usually omit the indices and write V for Vv,n, Q for Qy^n respectively. 
To simplify notations we will often also work with representatives instead of equivalence 
classes [$] G Q. 

The interpretation of the Grassmann manifold as equivalence classes of orthonormal bases 
spanning the same A^-dimensional subspace is just one way to define the Grassmann 
manifold. We can as well identify the subspaces with orthogonal projectors onto these 
spaces. To this end, let us for $ = ((/^i, . . . , yjAr) G denote by -D$ the L2-orthogonal 
projector onto span{(y9i, . . . , v^at}- It is straightforward to verify 

Lemma 1. There is a one to one relation identifying Q with the set of rank N L2- 
orthogonal projection operators . 

In the following, we will compute the tangent spaces of the manifolds defined above for 
later usage. 

Proposition 1. The tangent space of the Stiefel manifold at^EVis given by 

W = {xev''\ = - G M^^^} . 

The tangent space of the Grassmann manifold is 

r^^G = {Vr G V^l = G M^^^} 



(span{v5i,. . . ,V57v}" 



Thus, the operator (X — 1^$), where is the L2-projector onto the space spanned by $, 
is an L2- orthogonal projection from V'^ onto the tangent space T\^Q . 
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Proof. If we compute the Frechet derivative of the constraining condition 

^($) := ($^$) -1 = 

for the Stiefel manifold, the first result follows immediately. To prove the second result, 
we consider the quotient structure of the Grassmann manifold and decompose the tangent 
space 7$V of the Stiefel manifold at the representative $ into a component tangent to 
the set [$], which we call the vertical space, and a component containing the elements of 
7$V that are orthogonal to the vertical space, the so-called horizontal space. If we move 
on a curve in the Stiefel manifold with direction in the vertical space, we do not leave 
the equivalence class [$]. Thus only the horizontal space defines the tangent space of the 
quotient Q = V/0{N). The horizontal space is computed in the following lemma, from 
which the claim follows. □ 

Lemma 2. The vertical space at a point $ G V (introduced in the proof of proposition^ 
is the set 

{<I>M|M = -M^ G M^^^}. 
The horizontal space is given by 

{W G {W^^) = G M^^^} . 

Proof. To compute the tangent vectors of the set [$], we consider a curve c{t) in [$] 
emanating from $. Then c is of the form c(t) = $U(t) for a curve U(t) G 0{N) with 
U(0) = Itvxtv- Differentiating I^xiv = U{t)U{tf at t = yields U'(0) = -U'(0)^ and we 
get that every vector of the vertical space is of the form $M where M is skew symmetric. 
Reversely, for any skew symmetric matrix M we find a curve U(t) in 0{N) emanating 
from $ with direction M, and c(t) := $U()f:) is a curve with direction c(0) = $M, and 
thus the first assertion follows. 

To compute the horizontal space, we decompose W G 7$V into W = $M + W±, where 
W± ■.= W - ^ {^'^W) G M := (<I>^W^). Then M is an antisymmetric matrix, which 
implies that $M is in the vertical space, and that the horizontal space is given by all 
{W± = ly - $ ($^iy) \W G %V}. Let us note that this set is the range of the operator 
(/ — V^). This operator is continuous and of finite codimension. If W± = W — ^ ($-^1^) 
is in the horizontal space, then 

Reversely, if W E with (W^^) = 0, then W is in %V and from (T - V^)W = 
W — ($^iy) = W we get that W is in the range of / — being the L2-orthogonal 
projection from onto the tangent space T[$]^. □ 

To end this section, let us prove a geometric result needed later. 
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Lemma 3. Let [^] G Q, D* the L2-projector on span[\[^], V* is its expansion as above 
and \\.\\ is the norm induced by the L2 inner product. For any $ = G V 

sufficiently close to [\&] & Q in the sense that for all i E {!,. . . ,N}, ||(/ — D*)(pi\\ < 6, 
there exists an orthonormal basis ^ &V o/span[\l/] for which 

Proof. For i = 1, . . . , N , let 

ijji = aicgmm{\\ip-Lpi\\,ip espan{ipi\i = l,...,N},\\ip\\ = l} = D*(^i/\\D*Lpi\\, 

and set \1/ := {ipi, . . . ,iPn)- If we denote by Pi the L2 projector on the space spanned by 
it is straightforward to see from the series expansion of the cosine that 

{I-D*)ip, = {I-Pi)^^ = ifi-ip^ + 0{\\{I - D*)ip,,\\^) (2.3) 

The fact that \E' ^ V is remedied by orthonormahzation of \l/ by the Gram-Schmidt 
procedure. For the inner products occurring in the orthogonahzation process (for which 
i 7^ j), there holds 

= - {{I - D*)^,,i^,) - {{I - D*)^„ (I - D*)^j) + Om - D*)^,\\'). 
= 0(\\(I-V*)M^) 



where we have twice replaced (fi — ipi by (/ — D*)(pi according to (2.3) and made use 
of the orthogonality of D*. In particular, for $ sufficiently close to [\E'], the Gramian 
matrix is non-singular because the diagonal elements converge quadratically to one while 
the off-diagonal elements converge quadratically to zero. By an easy induction for the 
orthogonahzation process and a Taylor expansion for the normalization process, we obtain 
that \E' differs from the orthonormalized set ^ := {ipi, . . . only by a error term 

depending on ||(/ — X'*)$|p. Therefore, 

^r-A = + Om-V*M'') = (I-D*)^, + om~v*)^\\'), 

so that 

= (/-!?*)$ + 0(||(/-I?*)$||2), 
and the result is proven. □ 
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2.3 Optimality conditions on the Stiefel manifold 



By the first order optimality condition for minimization tasks, a minimizer [\E'] E Q oi the 
functional : ^ — M, $ i— > J7'($) over the Grassmann manifold Q satisfies 

= for all 5$ G , (2.4) 

i.e. the gradient J7''(\E') G (V')^ = {V'^)' vanishes on the tangent space Tq,Q of the 
Grassmann manifold. This property can also be formulated by 

{{S^fj'i^)) = for all (5$ G %]^, 

or equivalently, by Lemma [T| 

(((X - V^)J'{'^), $)) = for all $ G , (2.5) 

that is, in strong formulation, 

{I-V^)J'{^) = J'{^)-^A = Oe{VY, (2.6) 

where A = {{iJ'i^))j,ipi))f^j=i and iJ'i^))i G V is the i-th component of J'i"^). Note 
that this corresponds to one of the optimality conditions for the Lagrangian yielded from 
the common approach of the Euler-Lagrange minimization formalism: Introducing the 
Lagrangian 

/:($,A) := ^ (j($) + ^A,,((^„^,)^^ -5„)) , (2.7) 

the condition for the derivative restricted to V^, here denoted by C^^''^\'^,A) for conve- 
nience, is given by 

TV 

(VI/, A) = j'(^q;)-(J2Kk^Pk)l, = 0e{V'f. (2.8) 

k=l 

Testing this equation with ipjjj = Ij • • • > ^) verifies the Lagrange multipliers indeed agree 



with the A defined above, so that (2.5) and (2.8) are equivalent. Note also that the 
remaining optimality conditions, 

dC 1 

of the Lagrange formalism are now incorporated in the framework of the Stiefel manifold. 



From the representation (2.6), it follows that the Hessian £'^^'*^(\1/, A) of the Lagrangian 
(2.8), taken at the minimum \1/ and with the derivatives taken with respect to is given 

by 

£(2-*)(^^A)$ = J^"(^)$-$A. 
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As a necessary second order condition for a minimum, A)^^'*) has to be positive 
semidefinite on T[q,]Q. For our convergence analysis, we will have to impose the stronger 
condition on £^^'*)(\E', A) being elliptic on the tangent space, i.e. 

((£(2'*)(^,A)<5$ , 5$)) > 7 ||5$||yiv, for all 5$ G %]^. (2.9) 

It is an unsolved problem if this condition holds in general for the minimization problems 



of the kind (1.1) or if it depends on the functional under consideration; in particular, it 
is not clear whether it holds for the functionals of Hartree-Fock and density functional 
theory. In the case of Hartree-Fock, it suffices to demand that £'^^'*^(\Ef, A) > on Tj^j^ 
because this already implies £'^^'*)(\E', A) is bounded away from zero, cf [35J. For the 
simplified problem, we will show in Lemma |4] that the assumption holds for symmetric 
operators A fulfilling a certain gap condition. 



3 Minimization tasks in electronic structure 
calculations 

We will now particularize the results of the last section to the functionals common in 
electronic structure calculation. As the following section will show, the applications of 
interest in electronic structure calculations deal with the minimization of functionals 
for which the gradient can be written as J7''($) = Aq>^, where ^4$ : V ^ V (and A<j, its 



extention to V by (2.1)). We conjecture that if the functional J' only depends on the 



electronic density, that is, if condition (1.2) holds, this form of i7($) is valid in general, 
i.e. for each $ G ^, there is an operator ^4$ so that i7'($) = ^<i,$. Nevertheless, we 
decided to formulate the algorithms (except algorithm 3) for i7'($) rather than for Ai^, 
to emphasize the minimization viewpoint we pursue in this work and to display that the 
concrete structure of the Fock or Kohn-Sham operators does not enter anywhere in the 
proof of convergence given in section |5| 

In this section, we will remind the reader of some basic facts about Hartree-Fock and 
Kohn-Sham theory, where our emphasis will be on the ansatzes leading to the problem 



of minimizing a nonlinear functional (1.1). Also, we will review the concrete form the 



operator J7'($) = in (1.4) has in these applications. For a more detailed introduction 



to electronic structure calculations, we refer the reader to the standard literature [HI |T2l 



|26lll3]. At the end of this section, we will investigate the simplified problem (1.5) and its 
connection to eigenvalue computations. 
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3.1 Hartree-Fock and Kohn-Sham energy functionals 
in quantum chemistry 



The commonly accepted model to describe atoms and molecules is by means of the 
Schrodinger equation, which is in good agreement with experiments as long as the en- 
ergies remain on a level at which relativistic effects can be neglected. We are mainly 
interested in the stationary ground state of quantum mechanical systems, given by the 
eigenfunction belonging to the lowest eigenvalue of the Hamiltonian H of the system. In 
the Born-Oppenheimer approximation the Hamiltonian of the (time-independent) elec- 
tronic Schrodinger equation = E'if is given by 

^ N* N* M ^ N* ^ 

^ •~ r, 5^ 5^ 5^ |U. _ P II 9 Wr- — T • II 

1 = 1 1 = 1 U = l lJ = l^l^J ■> 

Here, A^* denotes the number of electrons, M the number of the nuclei, and Zy^ the 
charge respectively the coordinates of the nuclei, which are the only fixed input parameters 
of the system. Note that we use atomic units, so that no physical constants appear in the 
Schrodinger equation. We also neglect the interaction energy between the nuclei, since 
for a given constellation (i?i, . . . , -Rm) of the M nuclei this only adds a constant to the 
energy eigenvalues. Due to the Pauli principle for fermions, the wave function is required 
to be antisymmetric with respect to permutation of particle coordinates. It is easy to see 
that every such antisymmetric solution can be represented by a convergent sum of Slater 
determinants of the form 

V'lL(a;i,Si...,XAr*,SAr.) := ^l=det(v9i(a;j,Sj)), G M^ Si = ±^ 

where $ = (<^i)il*i G R^i^ x {±|})^* and = 6ij. In Hartree-Fock (HF) theory, 

one approximates the ground state of the system by minimizing the Hartree-Fock energy 
functional $ Jhf{^) '■= {HipsL,ipsi) over the set of all wave functions consisting of 
one single Slater determinant ipsL^^T-^ Si . . . ,xn*, sjy*). Additional simplification is made 
by the Closed Shell Restricted Hartree-Fock model (RHF), given in a spin-free formulation 
for = N*/2 pairs of electrons, so that $ = i<fi)fLi G H'^{M.^)^ =: V^. Abbreviating 
V{x) := — J2tLi II \ II , the corresponding functional reads 



1^=1 \\x-R^\\ ' 

N „ ^ ^ N 



2 



A minimizer of Jhf is named Hartree-Fock ground state. Its existence has been proven 
in the case that Y^^^i Z^> N -1 ([^H], 
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The energy functional of the Kohn-Sham (KS) model can be derived from the Hartree- 
Fock energy functional by two modifications: First of all, as a consequence of the Hohenberg- 
Kohn theorem (cf. [21]), it is formulated in terms of the electron density n{x) = \ P 
rather than in terms of the single particle functions; secondly, it replaces the nonlocal and 
therefore computationally costly exchange term in the Hartree-Fock functional (i.e. the 



fourth term in (3.1)) by an additional (a priori unknown) exchange correlation energy 
term Exc{n) also depending only on the electron density. The resulting energy functional 
reads 



N 



i=l Tnio mr, r^o ^ 



Determining the ground state energy of Kohn-Sham theory then consists in a minimiza- 
tion of J'ks over all $ = {ipi, . . . ,ipN) G V'^ with {ipi,ipj) = 6ij. Since the exchange 
correlation energy E^c is not known explicitly, further approximations are necessary. The 
most simple approximation for E^c is the local density approximation (LDA, cf. [^13]) de- 
fined as E^^^{n) = J n{x)e^^^{n{x)) dx, where e^f"^ denotes the exchange- correlation 

energy of a particle in an electron gas with density n. If we split this expression in an 
exchange and a correlation part, we get 

E^^^^in) = E^^^\n)+E^^\n) = j n{x)e'^f\n{x)) dx + j n{x)e';:^\n{x))dx, 

RS R3 

where in the exchange part, e^^^iji) = —Con^ and Cd '■= f (f )^^^ is the Dirac constant. 
For the correlation part E^^^{n), the expression e^^^ is analytically unknown, but can 
be calibrated e.g. by Monte-Carlo methods. We note that a combination of both HF and 
density functional models, namely the hybrid B3LYP, is experienced to provide the best 
results in benchmark computations. 



3.2 Canonical Hartree-Fock and Kohn-Sham equations 

For the HF and KS functionals, we can compute the derivative of and the Lagrange 
multipliers at a minimizer explicitly. 

Proposition 2. For the functional J'hf of Hartree-Fock, Jhf{^) — -A-qt^ G iV')^ , where 
A<i, = F^^ : H^{E?) — > H~^{E}) is the so-called Fock operator and A<i> is defined by 



through ( 2. 1 ); using the notation of the density matrix 



p>s>{x,y) := N J ips^{x,X2, . . . ,xn) ^^^(y, X2, . . . , xa?) dx2---dxN 

R3(iv-i) 
N 

= '^Vi{x)vi{y) 



i=l 
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and the electron density n$(x) := p$(x,x) already introduced above. It is given by 

+ f^^dy^ix)- [ 'f^y^^!;y^ dy. 

J \\x-y\\ J \\x-y\\ 

For the gradient of the Kohn-Sham functional Jks, there holds the following: Assuming 
that Exc in Jks is differentiate and denoting by v^c the derivation of E^c with respect 
to the density n, we have J7''($) = .4$$ G {V')'^ , with A,^ = F^^ the Kohn-Sham 
Hamiltonian, given by 

F^^ifi := -^A(fi + V{x)ifi + (n'k^j(fi + Vxc{n)ipi. 



In both cases, the Lagrange multiplier A of (2.8) at a minimizer = {tjji, . . . , ^pN) is given 
by 

Xij = {Aii,ilJi,iljj). (3.2) 

There exists a unitary transformation U G 0{N) amongst the functions ipi, i = 1, . . . , N 
such that the Lagrange multiplier is diagonal for \[^U = {ipi, . . . , iPn), 



Xij := {Aipi,ipj) = Xi6, 



so that the ground state of the HS resp. KS functional (i.e. minimizer of J') satisfies the 
nonlinear Hartree-Fock resp. Kohn-Sham eigenvalue equations 



F^''4ji = X,^lji, resp. F^^^^ = Xi'^j,, A, G M, i = l,...,N, (3.3) 

for some Ai, . . . , A^v G M and a corresponding set of orthonormalized functions \1/ = {ipi)f=i 
up to a unitary transformation U. 

The converse result, i.e. if for a collection $ = ((/^i, . . . ^^Pn) belonging to the lowest 



eigenvalues of the Fock operator in (3.3), the corresponding Slater determinant actually 
gives the Hartree-Fock energy by J7'($) = {Hipgi^, ipsi)^ is not known yet. 



3.3 Simplified problem 



The practical significance of the simplified problem (1.5) is given by the following result, 
which shows that for symmetric A, the minimization of J'a is indeed equivalent to finding 
an orthonormal basis {ipi : 1 < i < N} spanning the invariant subspace of A given by the 
first eigenfunctions of A. 
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Proposition 3. Let A in the simplified problem (1.5) a bounded symmetric operator. The 
gradient of the functional J a is then given by J7''($) = G iV')^ . Therefore, is a 
stationary point of L if and only if there exists an orthogonal transformation U such that 
'^\J = {ipi, . . . jIPn) G consists of N pairwise orthonormal eigenf unctions of A, i.e. 



Aipk = ^ki^k for k = 1, . . . , A^; in this case, there holds J^i^^) = Ylk=i ^k- The 
of J is attained if and only if the corresponding eigenvalues Xk, k = 1, . . . , N are the N 
lowest eigenvalues. This minimum is unique up to orthogonal transformations if there is 
a gap Xn+i — Xn > 0, so that in this case, the minimizers \1/ = argmin J are exactly 
the bases of the unique invariant subspace spanned by the eigenvectors according to the N 
lowest eigenvalues. 



minimum 



3.4 Comparison of direct minimization and 
self consistent iteration 

Self consistent iteration consists of fixing the Fock operator F*^") = F^(n) for each iterate 
the simplified problem is then solved in an inner iteration loo for A = F^^^; the solu- 
tion $ defines the next iterate $("'+^) of the outer iteration, by which the Fock operator is 
then updated to form F^"'~^^\ defining the simplified problem for the next iteration step. 
For the solution of the inner problems with a fixed Fock operator, Proposition |3] from the 
last section applies and the algorithms presented in the next section can be used. Self 
consistent iteration is faced with convergence problems though, which can be remedied 
by advanced techniques: With an appropriate choice of the update, the ODA-optimal 
damping algorithm |^^, convergence can be guaranteed. 



Direct minimization corresponds to the treatment of the nonlinear problem (1.1) for the 
Hartree-Fock or Kohn-Sham functional with the gradient algorithm 1 from the next sec- 
tion. Direct minimization thus differs from the self consistent iteration only in that the 
Fock operator is updated after each inner iteration step. Therefore, direct minimization 
is preferable if the update of the Fock operator is sufficiently cheap. This is mostly the 
case for Gaussians but not for the plane wave or wavelet basis or finite differences. 



4 Algorithms for minimization 

In this section we will introduce three related algorithms to tackle the minimization 



problem (1.1) in a rather general form. Their convergence properties will be analysed 



in the next section. 
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4.1 Gradient and projected gradient algorithm 



We will consider a gradient algorithm for the constrained minimization problem; the 
motivation for this is given by the following related formulation (cf. [32] for this concept): 
With an initial guess ^^^^ E V, [$*^°^] G the gradient flow on V, resp. Q is given by the 
differential 

- j'im), m = \/6<^e 7j^(,)]^. (4.1) 

Using the fact that X — is projecting onto the tangent space T^^-^Q, this algebraic 
differential initial value problem can be rewritten by an ordinary initial value problem for 
the gradient flow on V, 

or, equivalently, 

where the bracket [., ..] denotes the usual commutator. Denoting by (J7''($(t)))i the i-th 
component of the gradient J'{^(t)) and letting A(t) := [{{J'{^{t)))i,(pj{t)))^^^^, we 
obtain the identification 

which we will make use of later. 

There holds ^^^^ ^ for t ^ oo, so we are looking for the fixed point of this flow 



\1/ = limt^oo ^{t) rather than its trajectory. Equation (4.3) suggests the projected gradient 



type algorithms presented below. In algorithm 1, corresponding to an Euler procedure for 



the differential equation (4.1), the gradient at a certain point ^{t) is kept fixed (and being 
preconditioned) for non-differential stepsize, so that the manifold is left in each iteration 
step. Therefore, a projection on the admitted set is performed in each iteration step. 
Note also that the role of the preconditioners is crucial, see the remarks following 
algorithm 1. 

Algorithm 1: Projected Gradient Descent 

Require: Initial iterate <i>^°) G V ; 

evaluation of J^' {^^^^^) and of preconditioner(s) B^^ (see comments below) 
Iteration: 
for n = 0, 1, . . . do 

(1) Update A(") := ( J^'($(")), e M^^^, 

(2) Let $("+1) - B-\j'{^^"'>) ~ $(")A(")), 

( = - B-\A^(^)^^"'> - $(")A(")) for the case that J'{^) = A^^.) 

(3) Let = p$("+i) by projection P onto V resp. G 
endfor 
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Some remarks about this algorithm are in order. First of all, note that if Algorithm 1 is 
applied to the ansatzes in electronic structure calculation as portrayed in section |3| the 
gradient is given by = with the Fock- or Kohn-Sham operator or 

a fixed operator Aq, = A for the simplified problem. Therefore, {J''(^^"'^ — $(")A*^"^)j = 
y4<j,(n)0j-"^ — X]jLi(^<i>(")0i"''5 is the usual "subspace residual" of the iterate 

which is a crucial fact for capping the complexity of the algorithm in section [6j 



Next, let us specify the role of the preconditioner B^^^ used in each step. This precondi- 



tioner is induced (according to (2.1)) by an elliptic symmetric operator Bn ■ V ^ V, 



which we require to be equivalent to the norm on in the sense that 



~ Ml, yipeV = H\R'). (4.5) 

For example, one can use approximations of the shifted Laplacian, B ^ a(— lA + C), 
as is done in the BigDFT project. This is also a suitable choice when dealing with plane 
wave ansatz functions using advantages of FFT, or a multi-level preconditioner if one has 
finite differences, finite elements or multi-scale functions like wavelets |25l [3 |16[ [3]. 

For the simplified problem, the choice B^^ = aA^^ corresponds to a variant of simulta- 
neous inverse iteration. The choice 

corresponds to a simultaneous Jacobi-Davidson iteration. 

To guarantee convergence of the algorithm, the preconditioner B chosen according to 
the guidelines above also has to be properly scaled by a factor a > 0, cf. Lemma |6| 
The optimal choice of a is provided by minimizing the corresponding functional over 
span <|)('^+^)} (a line search over this space), which can be done for the simplified 

problem without much additional effort. For the Kohn-Sham energy functional, it will 
become prohibitively expensive. However, line search and subspace acceleration like DIIS 
[37] will improve the convergence speed. Note that in this context, one might as well use 
different step sizes for every entry, i.e. i3$ = {aiB{pi, . . . , a^BipN). 

Next, let us make a remark concerning the projection onto Q. It only has to satisfy 
span {^["''^^^ : 1 < < N} = span : 1 < < N}. For this purpose any orthogo- 

nalization of {^j-"^^'* : 1 < ^ < A^} is admissible. For example, three favorable possibilities 
which up to unitary transformations yield the same result are 

• Gram-Schmidt orthogonalization, 

• Diagonalization of the Gram matrix G = {{^i^^^K 0^f'^^^))i^j=i by Cholesky factor- 
ization, 
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• (For the problems of section |3| i.e. where J'{^) = A^^:) 

Diagonahsation of the matrix A${„+i) := {{Ai^(n)ipl"'^^\0^^~^^^))f^j^i by solving an 
N X N eigenvalue problem. 

Parallel to the above algorithm, we consider the following variant in which the descent 
direction is projected onto the tangent space 'Zj-(j,(n)j Q in every iteration step. It will play an 
important theorectical role considering convergence of the local exponential parametriza- 
tion, i.e. algorithm 3. 

Algorithm 2: Modified Projected Gradient Descent 

Require: see Algorithm 1 
Iteration: 

for n = 0, 1, . . . do 

(1) Update A(") := ( e K^>^^, 

(2) Let -(X-X'<j,(„))S-i(j^'($("))-$(")A(")), 

( = $(") - 5-H-^<i.(")*^"^ - $(«)A(")) for the case that J'{<^) = A^^) 

(3) Let = p$(«+i) hy projection P onto V resp. G, 
endfor 

Note again that the algorithms are given in a general form, where the preconditioner (or 
the corresponding parameter a, e.g. obtained by a kind of line search) may be chosen in 
each iteration step. In our analysis, we will consider a fixed preconditioner Bn = B m 
every iteration step, for which we will show linear convergence without further line search 
invoked. Thus, our analysis is in a way a worst case analysis for the algorithms under 
consideration. See also section [6] for improvements on the speed of convergence. 

4.2 Exponential parametrization 

Instead of projecting the iterate onto the Grassmann manifold Q in every iteration 
step, we will now develop an algorithm in which the iterates remain on the manifold with- 
out further projection. This will be achieved by following geodesic paths on the manifold 
instead of straight lines in Euclidean space, which has the advantage that during our cal- 
culations we do not leave the constraining set at any time so that no orthonormalization 
process is required. To apply the result of proposition |4| we will for this algorithm limit 
our treatment to the case where J7''($) = is given by a linear operator (see the 

discussion after algorithm 1). 

Recall that a geodesic is a curve c on a manifold with vanishing second covariant deriva- 
tive, i.e. 

j^c{t) := 7r,(t)c(t) = foralH, (4.6) 
where Tic(t) denotes the projection onto the tangent space at the point c{t). 
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Proposition 4. For any operator X : V ^ V for which G T\q,-\Q (where as always, 



X is defined by X by (2.1)), the antisymmetric operator 



X = {I - D^)XD^- D^X\I - D^), (4.7) 

satifies ^Y<I> = and c(t) := exp(t^)$ is a geodesic in Q emanating from point $ with 
direction c(0) = X^. 

Proof. The proof is straightforward; application of the projection equation yields (^c(t)) = 
(X-P,(i))c(t) = 0. □ 

If we now let, for any iterate 

= (/ - D^^^))B-\I - (4.8) 

the curve 

c{t) := exp(-t;^ ("))$(") 



with A"*^") from (4.7) is by the previous Lemma a geodesic in Q with direction 

c(0) = -(/-D$(„))5-i(J-D^(„))A$(„)<l>(") 

which equals the (preconditioned) descent direction of the projected gradient descent 
algorithm of the preceding section. If we now choose the next iterate as a point on this 
geodesic, we get the following algorithm: 

Algorithm 3: Preconditioned exponential parametrization 

Require: see Algorithm 1 
Iteration: 

for n = 0, 1, . . . do 

\> Follow a geodesic path on the Grassmann manifold with stepsize a, 



:= exp(-ai'("))$(")(with X from (|£8| and X defined by X via (|4J|) 

endjor 



Note that a similar algorithm. Conjugate Gradient on the Grassman Manifold, has al- 
ready been introduced in page 327. That paper also included numerical tests for a 
model system. The algorithm was also tested for electronic structure applications very 
different from those of the BigDFT program in [38]. A similar approach using the density 
matrix representation for electronic structure problems was also proposed in |12], where 
the authors move along the geodesies in a gradient resp. Newton method direction with- 
out preconditioning. 
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Like in this work, the stepsize a may be calculated in each iteration step using line search 
algorithms like backtracking linesearch or quadratic approximations to the energy term 
|36j . These often time consuming line searches may be omitted though if we choose a 
suitable preconditioner B = and set the stepsize a = 1 once and for all. 
The efficiency of this algorithm strongly depends on the computation of matrix exponen- 
tials needed to follow geodesic paths on the Grassmann manifold. A variety of methods 
can be found in [33], see also [H] for an analysis of selected methods. For some of these 
algorithms, there exist powerful implementations like the software package Expokit 
which contain both Matlab and Fortran code thus supplying a convenient tool for numer- 
ical experiments. 



5 Convergence results 

5.1 Assumptions, error measures and main result 

In this section, we will show linear convergence of the algorithms of the last section under 
the ellipticity assumption [T] given below. Additional results we give include the equiva- 
lence of the error of $, measured in a norm on V , and the error of the gradient residual 
(X — D)jr'(<l>), and quadratic reduction of the energy error 

Recall that in our framework introduced in section [2| we kept the freedom of choice to 
either use V := = H^{M.^), equipped with an inner product equivalent to the inner 
product (., for analysing the original equations, or to use V = Vh C as a finite 
dimensional subspace for a corresponding Galerkin discretization of these equations. In 
practice, our iteration scheme is only apphed to the discretized equations. However, the 
convergence estimates obtained will be uniform with respect to the discretization param- 
eters. The main ingredient our analysis is based on is the following condition imposed on 
the functional J', cf. section 2.3[ 



Assumption 1. Let \E' a minimizer of (1.1). The Hessian £'^^'*)(\E', A) : {V')^ of 



the Lagrangian £(\E',A) (given by (2.8)), where the derivatives are taken with respect to 
\E', is assumed to he -elliptic on the tangent space, i.e. there is 'y > so that 

((£(2''^)(^,A)5$ , 5$)) > 7 for all 5$ G 7|*]^. (5.1) 

Note again that for Hartree-Fock calculations, verification of £'^^'*)(\E', A) > on Tj^j^ 
already implies a), cf [35]. 



From section 
and only if 



2.2 



we recall that a)<I) = J"{'^)^ - <I>A, so that (5.1) is verified if 



((jr"(^)5$ - 5$A , (5$)) > 7 ||(5$||yiv, for all 5$ G (5.2) 
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holds, where A = (((i7'(^))j, V'j))i^=i above. From the present state of Hartree-Fock 
theory, it is not possible to decide whether this condition is true in general; the same 
applies to DFT theory. For the simplified problem, the condition holds if the operator A 
fulfils the conditions of the following lemma. 

Lemma 4. Let A : V V, ip t— > Aip a hounded symmetric operator, such that A has N 
lowest eigenvalues \i < ■ ■ ■ < \n satisfying the gap condition 

\n < inf{A I A e cx{A)\{X,, A^v}}- (5.3) 

Then assumption^ holds for the simplified problem (1.5). 



Proof. We estimate the two terms of (5.2) separately. Let us denote A = inf{A | A G 
cr(A)\{Ai, . . . , Aat}}. To the first term, the Courant-Fisher theorem ([39]) applies com- 
ponentwise to give the estimate {{A6^,6^)) > \\\6^\\yN- For the second, choosing 
U = {ui,j)fj^i G 0{N) so that U^AU = diag(Aj)^i, where Aj are the lowest N eigenval- 
ues of A, gives 

N N N 

((5$A,5$)) = (((5$(UU^AUU^),5$)) := 5^(^M,-,A,%,^Mfc,5(^,) 

i=l j=l k=l 

N 

= X] ^j^j,k{S'~Pj,Sipk) < AAr||5$||y]v. 
j,k=l 



SO that £(^)(\E', A) is elliptic on T^Q by the gap condition (5.3). □ 



To formulate our main convergence result, we now introduce a norm ||.||yjv on the space 
V^^, which will be equivalent to the (iy^)^-norm but more convenient for our proof of 
convergence. We will then state our convergence result in terms of these error measures. 

Lemma 5. Let B : V V the preconditioning mapping introduced in section^ so that 
in particular, B is symmetric and the spectral equivalence 

< {Bx,x) < Q\\x\\li 

holds for some < < Q and all x ^ V. Let us consider the mapping 

B-^:V'^V, B-^ := {I - D)B-\I - D) + D, (5.4) 

where D = projects onto the sought subspace. Then the inverse B satisfies {Bip, ip) = 
{ip, Bip) for all (p,ip E V , and for the induced B-norm \ \.\\^ on V there holds 

{B<p,<p) ~ ||<^||^i. 
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Using the notation (2.1), a norm on V'^ is now induced by the ||.||^-norm by 

ml^:={{m,^)) . (5.5) 

Note that this norm, as any norm defined on in the above fashion, is invariant under 
the orthogonal group of M^^^ in the sense that 

||$U||yJV = ||$||yiv (5.6) 

for all U G 0{N). In the Grassmann manifold, we measure the error between [$(i)], [$(2)] £ 
^ by a related metric d given by 

d{[^ii)],[^{2)]) ■■= ^mf^^||<l>(i)-<l'(2)U||y^. 



If [$(2)] is sufficiently close to [$(1)] G ^ it follows from Lemma [3] that this measure given 
by d is equivalent to the expression 

||(J-D^(J$(2)||v-, (5.7) 

in which we used the L2-orthogonal projector "P^^jj onto the subspace spanned by 
In the following, let us use the abbreviation D = for the projector on the sought 
subspace, whereever no confusion can arise. An equivalent error measure for the deviation 
of $ G V from the sought element \I' G V is then given by the expression 

||(J-P)<l>||v^.v, (5.8) 

which will be used in the sequel. In terms of this notation, our main convergence result 
is the following. 



Theorem 1. Under the ellipticity assumption (5.1), the following holds for any of the 
three algorithms formulated in section^ For ^'^^^ G Us{^) sufficiently close to ^, there 
is a constant X < 1 such that for all G Nq, 

||(J_l))$(n+l)||^^ < ;^.||(J-I?)$W||^^. (5.9) 

The rest of this section will be mainly dedicated to the proof of this theorem. For the sake 
of clarity, let us first sketch the proof to be performed: We will exploit the fact that the 
iteration mapping can be written in the form ^ - B-^{I - V^(-a))J'{¥''^) and 
is thus a perturbation of the mapping 

$W ^ $W _ B-^{I-V^)J'{<^^"'y). The estimate 
then splits in two main parts: The first will be a linear part incorporating the Hessian 
of the Lagrangian and the task will be to show that application of this linear part to an 



iterate <1>'^"'' G Q indeed reduces its error in the tangent space of \E' (as defined by (5.8)); 
here, our ellipticity assumption enters as main ingredient. The second part consists of 
showing that the remaining perturbation terms (including those resulting from projection 
on the manifold) are of higher order and thus asymptotically neglectable; the main lemmas 
entering are Lemma |3] above and Lemma [8] to be proven below. 
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5.2 Ellipticity on the tangent space 



In this section, we will first formulate a rather general result about how ellipticity on 
subspaces can be used to construct a contraction on these spaces and then specialize this 
to the tangent space at the solution \E' and assumption [T] in the subsequent corollary. 
Finally, we will then prove that our assumption 5A entering here is indeed true for the 
simplified problem (1.5). 

Lemma 6. Let W G G G W a Gelfand triple, U G W a closed subspace of W and 
S,T' : W ^ W two bounded elliptic operators, symmetric with respect to the G-inner 
product {.,..)g, satisfying 

^\\x\\w < {Sx,x)g < r||x||^, (5.10) 
and ^\\x\\w < {T'x,x)g < Q\\x\\w (5-11) 

for all X G U. Moreover, let S, T' both map the subspace U to itself. Then there exists a 
scaled variant T = aT' , where a > 0, and a constant (3 < 1 for which 

\\{I-T-'S)x\\t < P\\x\\t., (5.12) 
for all X G U, where := {Tx, x)g is the inner product induced by T. 

Proof It is easy to verify that for (3 := (16 -7??)/ (16 + 77?) < 1 and a := ^(1/^9 + 7/6) 
there holds 

-T-^S)x,x)t\ < P\\x\\l for all a; Gf/. (5.13) 



Due to the symmetry of T, S as mappings U ^ U, the result (5.12) follows. □ 



Let Xi,i = 1,...,N the lowest eigenvalues of A, ifji^i ■ 
eigenf unctions, and 

Vo = span {ipi-. i = 1,. . . ,N} 



1, . . . ,N, the corresponding 



(5.14) 



By Lemma [T| there holds (V^^)^ = '/[^j^, where = {ipi, . . . ,ipiy). The following corollary 
is the main result needed for estimation of the linear part of the iteration scheme. 



Corollary 1. Let J fulfil the ellipticity condition (5.1) and B' : V —> V a symmetric 



operator that fulfils (5.11) with T' = B' . Then there exists a scaled variant B = aB' , 



where a > 0, for which for any 5$ G T[^]Q there holds 

\\6<t>-B-\l-V)C'^^''^\'$,A)6<l>\\v^ < /5 



where (3 < 1 and B is defined by B via (5.4)- 
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Proof. Note that the restriction of B' is a symmetric operator Vq-^ —>■ V^^, so that the same 
holds for the extension B' as mapping T[q,]Q — > Tjij,]^. (X — V)C^'^''^^ also maps Vf^ — » 
symmetricly, so Lemma [6] applies. □ 



5.3 Residuals and projection on the manifold 

For the subsequent analysis, the following result will be useful. It also shows that the 
"residual" (T — r'^{„))j7''($*^"^) may be utilized for practical purposes to estimate the 
norm of the error (/ — D)^^"'\ 

Lemma 7. For 5 sufficiently small and ||(X — < 5, there are constants c, C > 
such that 

c||(j -!?)$(") 11^^ < ||(X-I?^,„))J'($("))||(v^), < C||(J-I?)<l>(")||yiv. (5.15) 
An analougeous result holds for gradient error — P)J''($^"^)||(yivy. 



Proof. Let us choose \E' G [^1^] according to Lemma [s] (applied to $ = $^")). Letting 
:= — ^, there holds by linearization and Lemma [s] (recall that we let D = D^) 

(X-P^m)J^'(<I>(")) = {I-V)J'{^) + (J-I?)£(2-i')(^^A)A^ + C(||(X-P)$(")||2 

= (j-i?)/:(2>*)(i|r^A)(j_p)$w + c'(||(x_p)$wii2 



yN 



w 



By assumption |T) ||(J-I?)/:(2,i')(^^ A)(J-I?)<l>(")||(y^y ~ 1 1 (J from which 
the assertion follows. The assertion for ||(X — T))J''(^^"'^)\\(^yny follows from the same 
reasoning by replacing £*^^'*)(\l/, A) by J^"{^) in the above. □ 

The last ingredient for our proof of convergence is following lemma which will imply that 
the projection following each application of the iteration mapping does not destroy the 
asymptotic linear convergence. 

Lemma 8. Let <|)("+^) = (0i,...,07v) the intermediate iterates as resulting from itera- 
tion step (2) in algorithm 1 or 2, respectively. For any orthonormal set $ G V fulfilling 
span[$] = span[<l>'^"+^^], its error deviates from that 0/$'^""'"-^^ only by quadratic error term: 

Wil-VMyM = ||(J-D)l>("+i)||y^ + 0(||(J-D)l'(")||^^) (5.16) 



Proof. First of all, note that if (5.16) holds for one orthonormal set $ with span[$] = 
span[$("'"'"^''], it holds for any other orthonormal set $ with span[$] = span[$*^""'"^)] because 
{I~V)^lJ\\vN = \\{I-V)^\\vN for all orthonormal U G 0{N). Therefore, we will show 



(5.16) for $ = (</?!, . . . , lpn) yielded from $("+^) by the Gram-Schmidt orthonormalization 
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procedure. Denote (^j = ipi'' + r-"-*, where for = B^^{X — X'^(„))j7''($''"^), we set 
Tj-"'* = s-""* or Tj-"^ = (/ — D${n))s-"'' for algorithm 1 or 2, respectively. From the previous 
lemma, we get in particular that < ||(-^ ~ for both cases (remember 

that D = Dq,). With the Gram-Schmidt procedure given by (p'j^ = (p^ — J2j<i{Vk,'^j)'^j, 
Vk = v'k/\Wk\\i lemma is now proven by verifying that in each of the inner products 
involved, there occurs at least one residual and that, on top of this, for the correc- 

tion directions there holds (/ - D)^'^ = 0(||(X - r')$(")||yiv) + C(Ei<fc = 
C>(||(T- P)$(")||yiv). Therefore, the correction terms are of 0(||(J - 'D)l>(")||^jv), thus 
proving ip'/^ — 0k = C'dK/ — D)^\\ypf). It is easy to verify that the normalization of ip'/^ 
only adds another quadratic term, so the result follows. □ 



5.4 Proof of Convergence 



To prove (|5^ for Algorithm 1, we define = ^-B-^{I-V^)J'{^), so that = 

p(jr($W)), where P is a projection on the Grassmann manifold for which [P(^($(")))] = 
For fixed n, let us choose ^ G span[\l/] according to Lemma [sj so that, using 
the abbreviation V := V^, 



< 



(J-P)<l>(") + C(||(J-P)<I> 



)||2, 



(n)||2 



Introducing A\E' := — \E', there follows by linearization 



Lemma [s] 



||(J_1)) $("+!) 11^^ 

II (J - I?)^(<l>(")) II V'iv + 0( II (X - f^^.) 

||(J-I?)J^(^') + (J-I?)r(^)A^||yiv +C(||(J-P)<I>(")||^ 

||(J 11^^ +0(||(J _!?)$(«) ||2^^) 

II (J -V){I- B-\l - V)C^^^^\^, A)) (J - I?)<l>(") ||y^ 
+0(||(J-D)$(")||2^.) 



(5.17) 
(5.18) 



(5.19) 
(5.20) 
(5.21) 
(5.22) 
(5.23) 



where we have used (5.18) and the fact that (X— P)jF(\E') is zero. The proof is now finished 
by noticing that 

(J -v){x-B-\x- r')/:^^'*) (^, A)) ( J - D)^ 
= {l-B-\l-V)C^^''^\^,k)YT-V)^, 

so that corollary [T] applies to give 

< ?9||(J-I?)<l'(")||yiv + C»(||(J-P)<l>(")||^^) < xll(2:-i^)$^"^| 
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where x < 1 for ||(X — I?)$*^"^||yiv small enough to neglect the quadratic term. □ 



The convergence estimate (5.9) for Algorithm 2 is easily derived from this: Consider 



J^2($) = $ - (J - V^)B-\X - (5.24) 

for which <|)("+^) = P(JF2(<I>("))) for the iterates of Algorithm 2. Differentiation of JF2 at ^ 
chosen as before gives 

= J- (J-I?)S-i(J-P)/:(2)(i|r^A)A^ + C(||(J-r')$(")||^^), 

(note that derivation of the projector on the left hand side with respect to 'I' results 
in a zero term), so that the same reasoning as above gives 



(X-I?)$("+^)| 



yN 



< \\(J-V){I- B-\l - P)/:(2'*)(^, A)) (X - V)^\\yN + 0(|| (X - !?)$(" 



)||2 



< xll(x-P)<f(")| 

with X < 1 fo^^ close enough to □ 

To prove the convergence of the exponential parametrisation (Algorithm 3) defined by 

$(n+i) ._ (^_aX^ 

it is enough to notice, cf. the remarks after Lemma |4| that we follow a geodesic path in 
direction (/ — 'r'^(n))i3~^(y4$(n)$*^"^ — ^("^A*^'^)), which is equal to the descent direction 
of Algorithm 2. Due to the definition of the tangent manifold, <|)("+^) again differs from 



JF2($'^"^) (defined by (5.24)) only by an asymptotically neglectable quadratic error term. 

□ 

5.5 Quadratic convergence of the energy 

For the Rayleigh quotient i?((/)*^")), i.e. for the simplified problem and = 1, it is known 
that -R(0''"'*) — -R(V') ^ IIV^ — ^'■"'■'lly- To end this section, we will show that this property 
holds also for the computed energies, provided that the constraints are satisfied exactly 
and the functional is sufficiently often differentiable. The latter is only known for Hartree- 
Fock and the simplified problem. Since the exchange correlation potential is not known 
exactly, this question remains open in general for the density functional theory. 
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Theorem 2. Provided that J is two times differentiable on a neighborhood Usi"^) C 
of the minimizer , and that for fixed ^ G Us{^), J" is continuous on + (1 — G 
[0,1]}, the error in the energy depends quadratically on the approximation error of the 
minimizer i.e. 

Ji<^)-jm < . (5.25) 

Proof. Let us choose a representant of the solution \I/ according to Lemma [3j Abbreviating 
e = $ - ^, we can use J^'(^)((J - !?)$) = to find that 

jr'(^)(e) = J\^){{I-V)<I>) + 0(||(J-P)<l>|n = 0(||(T-P)$||2) 

so that 



1 

Ji<^)-jm = j J'{^ + se){e)ds + ]^J'm{e) 







By integration by parts, 



1 1 
^(/(O) + /(!)) = j f{t)dt + lis-^-)ns)ds, 





so that 



i 

= l((y($),$-v[/)) _ J(s-^)J"{^ + se){e,e)ds + Oi\\{I-V)<^\ 



For estimation of the first term on the right hand side, recall from (5.15) that 

\\{i-v)j'm\y. < \\ii-v)ny., 

and therefore 

= o{\\{i~vm% 
1 

while for the second term, \ J{s - l)J"{^ + se){e,e)ds\ = 0{\\e\\^) = Oi\\{I - V)^\\^) 



follows from the continuity of Jj" and, again, the usage of Lemma [s] □ 
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6 Further Comments and Conclusions 



Before we conclude this article with numerical examples, we would like to make some 
comments about the complexity of the numerical schemes when applied to the problems 
of section [3} and about the potentialities for accelerating convergence of the iteration 
scheme. 

Complexity: Concerning disk storage, the task is to compute functions ip ^ ^h, so 
0{Ndim.Vh) memory is needed to store the orbital functions, while storage of the dis- 
cretization of the Fock operator A requires at most 0((dimV/j)^) in the general and worst 
case, but only 0{dimVh) for sparse discretizations. Regarding computational demands, the 
non-zero entries of a sparse discretization of A are of 0{dimVh), so that the complexity 
of the application of A depends linearily on dim\4. The computation of {A(/)[^~^^\ (f)^J^^^'') , 
and {(f)["''^^\(t>^J^'^^^) needs 0{N'^{dimVh)) operations in the case of sparse discretizations 
(and 0{N'^{dimVh)'^) in the worst case). The orthogonalization procedure, i.e. the pro- 
jection onto the Stiefel manifold usually has a complexity 0{N'^dimVh). To relate the 
above complexities to the size of the electronic system, it is also interesting to discuss 
how large dimVh^min has to be chosen for a given size N. To this end, we might fix a 
given maximal error e per atom or electron (usually requested to be smaller than the 
intrinsic modeling error of DFT or HP models) and determine the minimal ansatz space 
dimension dimVh,min{N) that keeps the numerical error under that error e. If we then 
consider the scaling of dim Vh,min with respect to the size of the system A^, it turns out that 
dimVh^rniniN) = 0{N), where the constant in front of is extremely large for systematic 
basis functions and surprisingly small for Gaussian type basis functions. Therefore, the 
natural scaling of the orbital based DFT and/or HP computations with respect to the 
size A^ of the underlying system gives an overall complexity of 0{N^) (or even 0{N^) for 
non-sparse discretizations). 

This can be improved if the discretization of the individual orbitals 0-"^ requires sub- 
stantially less than dimV/i DOFs. In an optimal case, one may archieve 0{1) for a fixed 
accuracy per atom; this is for example the case if the diameter of support of (f)['^^ is of 0{1), 
i.e. the support is local. In this case, the total complexity scales only linearly with respect 
to A^. Usually, the eigenfunctions xpi have global support. For insulating materials, though, 
there exists a representation \E'/oc such that [^loc] = [^] & Q and \ipioc,iix)\ < e~"'^~^'', 
a >> sufficiently large. These representations are called maximally localized or Wannier 
orbitals. Linear scaling 0{N) can be achieved if, during the iteration, the representant 
in the Grassmann manifold is selected and approximated in a way that the diameter 
of support is of C(l). This is the strategy pursued in Big DFT to achieve linear scaling, 
[T71 [22] . We defer the further details to a forthcoming paper. A related approach, com- 
puting localized orbitals in an alternative way was proposed by [6] and exhibits extremely 
impressing results. 
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Convergence and Acceleration: In the present paper we have considered linear conver- 
gence of a preconditioned gradient algorithm. For the simplified model, this convergence 
is guaranteed by the spectral gap condition, in physics referred as the HOMO-LUMO gap 
(i.e. highest occupied molecular orbital-lowest unoccupied molecular orbital gap). For the 
Hartree-Fock model, this condition is replaced by the coercivity condition The same 
condition applies to models in density functional theory, provided the Kohn-Sham energy 
functionals are sufficiently often differentiable. Let us mention that a verification of this 
conditions will answer important open problems in Hartree-Fock theory, like uniqueness 
etc. The performance of the algorithm may be improved by an optimal line search, re- 
placing B by an optimal Except for the simplified problem, where an optimal line 
search performed like in the Jacobi-Davidson algorithm as a particular simple subspace 
acceleration, optimal line search is rather expensive though and not used in practice. 
Since the present preconditioned steepest decent algorithm is gradient directed, a line 
search based on the Armijo rule will guarantee convergence in principle, even without a 
coercivity condition [5], [H] . 

In practice, convergence is improved by subspace acceleration techniques, storing iterates 
. . . , and compute $("+^) from an appropriately chosen linear combina- 

tion of them. Most prominent examples are the DIIS [3Z| and conjugate gradient [21 0] 
algorithm.The DIIS algorithm is implemented in the EU NEST project BigDFT, and 
frequently used in other quantum chemistry codes. Without going into detailed descrip- 
tions of those methods and further investigations, let us point out that the analysis in 
this paper provides the convergence of the worst case scenario. Second order methods, in 
particular Newton methods have been proposed in literature p5], but since these require 
the solution of a linear system of size NdimVh x NdimVh, they are to be avoided. 



7 Numerical examples 

The proposed direct minimization algorithm 1 is realized in the recent density functional 
code bigDFT [U], which is implemented in the open source ABINIT package, a com- 
mon project of the Universite Catholique de Louvain, Corning Incorporated, and other 
contributors jlHl [231 EH US]. It relies on an efficient Fast Fourier Transform algorithm 
[19] for the conversion of wavefunctions between real and reciprocal space, together with 
a DIIS subspace acceleration. We demonstrate the convergence for the simple molecule 
cinchonidine {CigH22N20) of moderate size A = 55 for a given geometry of the nuclei 



displayed in figure TA^. Despite the fact that the underlying assumptions in the present 
paper cannot be verified rigorously, the proposed convergence behavior is observed by 
all benchmark computations. The algorithm is experienced to be quite robust also if the 
HOMO-LUMO gap is relatively small. 
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Figure 7.1: Atomic geometry and electronic structure of cinchonidine 




Figure 7.2: Convergence history for the direct minimization scheme (left) and with DIIS accel- 
eration (right) for different mesh sizes. 



memory runtime 




Figure 7.3: Memory requirements (left) and computing time (right) for direct minimization 
algorithm with and without DIIS acceleration. 
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For our computations, we have used a simple LDA (local density approximation) model 
proposed by [20] and norm-conserving non-local pseudopotentials (21]. The orbital func- 
tions ipi are approximated by Daubechies orthogonal wavelets with 8 vanishing moments 
based on an approximate Galerkin discretization [18]. For updating the nonlinear poten- 
tial, the electron density is approximated by interpolating scaling functions (of order 16). 
The discretization error can be controlled by an underlying basic mesh size hgrid- 



In figure 7^, we demonstrate the convergence of the present algorithm for 4 different 
choices of mesh sizes, where the error is given in the energy norm of the discrete functions. 
The initial guess for the orbitals is given by the atomic solutions. Except in case of non- 
sufficient resolution {hgrid = 0.7), where we obtain a completely wrong result, convergence 
is observed. If the discretisation is sufficiently good, we do not observe much difference 
in the convergence history for different mesh sizes. Since the convergence speed depends 
on the actual solution, it is only possible to observe that the convergence is bounded by 
a linear rate. 

The number of iterations is relatively moderate bearing in mind that one iteration step 
only requires matrix-vector multiplications with the Fock operator and not a correspond- 
ing solution of linear equations. The DIIS implemented in BigDFT accelerates the iter- 
ation by almost halving the number of iterations and the total computing time at the 
expense of additional storage capacities, see also figure |7.2[ Further benchmark computa- 



tions have already been performed and will be reported in different publications by the 
groups involved in the implementation of BigDFT. 
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